On the Tail of the Overlap Probability 
Distribution in the Sherrington-Kirkpatrick 

Model 



Alain Billoire 

Service de Physique Theorique, CEA/DSM/SPhT 

Laboratoire associe au CNRS 
CEA-Saclay, 91191 Gif-sur-Yvette cedex (France) 

Silvio Franz 

ICTP, Strada Costiera 11, PO Box 563 
34100 Trieste (Italy) 

Enzo Marinari 

Dipartimento di Fisica, SMC and UdRl of INFM and INFN 
Universita di Roma La Sapienza 
P. A. Moro 2, 00185 Roma (Italy) 

billoirQspht . saclay .cea.fr 
f ranz@ictp . tr ieste . it enzo . marinariOromal . inf n . it 

February 1, 2008 

Abstract 

We investigate the large deviation behavior of the overlap proba- 
bility density in the Sherrington-Kirkpatrick model using the coupled 
replica scheme, and compare with the results of a large scale numer- 
ical simulation. In the spin glass phase we show that generically, for 
any model with continuous RSB, 1/JV log Pn(q) ~ — A {{\q\ — Qea)^ ■> 
and compute the first correction to the expansion of A in powers of 
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T c — T for the SK model. We also study the paramagnetic phase, 
where results are obtained in the replica symmetric scheme that do 
not involve an expansion in powers of q — q E A ° r T c — T. We give 
finally precise semi-analytical estimates of P(\q\ = 1). The overall 
agreement between the various points of view is very satisfactory. 

1 Introduction 

Although the solution of the Sherrington-Kirkpatrick (SK) model was 
proposed more than two decades ago [Q , and it physical interpretation 
is (now) quite clear |2j, some fundamental aspects of the model are still 
poorly understood. We are thinking here about the problem of the 
time scales of the model j3] , and to the question of chaos with respect 
to temperature variations [U[3], J us t to quote some recent works. 

In this article, we will revisit the question of the large deviation be- 
havior of the order parameter probability distribution P/v (q) of the SK 
model 013 El; using new analytical results, much improved numerical 
data and refinements in the data analysis. 

In the usual mean field approximation the probability distribution 
of the order parameter P/v(m) is read from the free energy consid- 
ered as a function of the order parameter. This method cannot be 
pursued in the replica approach, since the replica free energy has no 
physical interpretation for values of the order parameter different from 
the saddle point |2j. The coupled replica method was proposed in [g] 
to circumvent this problem and used in In short, the replica 

method is applied to a system of two identical copies of the original 
system with a constrained value of the mutual overlap q, with the 
result: 

^ log P N (\q\ > q EA ) = -A (\q\ - q E A? + O ({\q\ - q EA ) A ) , 

where q E A is the Edwards-Anderson order parameter (the maxi- 
mum allowed value for the overlap in the infinite volume limit). The 
proportionality coefficient A was computed in [§] to zeroth order in 
r = (T c 2 — T 2 )/2T C . In order to obtain a value of the coefficient A to 
compare with the results of our large scale simulation using the par- 
allel tempering algorithm, we compute here the first order term in r. 
We also present semi-analytical results for P(\q\ = 1). For T > T c we 
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present the results of an exact calculation, made in the replica sym- 
metric scheme, that is limited neither to small g's nor to the vicinity 
of T c . This calculation corroborates nicely the other analytical and 
numerical results. 

There are somehow two aspects in the comparison made between 
the analytical results and the Monte Carlo data. The first is to (cross) 
check both results. This is important since on one hand the analytical 
computation is very delicate and involves some unproven assumptions. 
The simulation is, on the other hand, limited to small systems finite 
number of disorder samples and finite statistics. Both approach are 
very complementary. The second aspect is in checking how much of the 
infinite volume physics (at low temperatures, in the spin glass phase) 
is already encoded in the sizes we can effectively study on our current 
state of the art computers (up to N = 4096 in the present case): we 
will see that, as far as the present computation is concerned, things 
do work well, and we are moving on solid ground. 



2 Large Deviations in Pn{q) f° r \q\ > 
Qea 



In this section we sketch the main steps and give the results of our 
computations of the function Pn(q) in the large deviation regime, i.e. 
for \q\ > qEA, for both T < T c and T > T c . 

The overlap probability distribution Pjv,j(9) depends on the dis- 
order sample J, and the function Pn(q) is obtained as Pn,j(q), where 
the over-bar denotes the disorder average. It is well known that in the 
thermodynamic limit the disorder average of Pn,j{q) is non- vanishing 2 
in the interval [0, ^sa]? while sample to sample fluctuations remain 
strong even for large systems The probability of the complemen- 
tary interval, i.e. of events with q > qEA, is exponentially small in 
the system size [Sj. For typical samples one has, independently of the 
disorder realization, 



where (3 = 1/T and F(q) is the self- averaging free energy cost of 

In the following we will take, without loss of generality, q > (assuming for example 
the presence of a suitable infinitesimal magnetic field). 

2 We assume that our starting Hamiltonian is symmetric under global spin reversal. 



PN,j(q) « P N (q) oc exp (—(3NF(q)) 
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keeping two replicas at a mutual overlap q: 



F(q) = —log ± £ e-«^W+^W) (ff - ^ E <™) • ( 2 ) 

It is quite natural to expect that in the thermodynamic limit P(q) = 
P(q). In an alternative scenario the fluctuations of Pjv j(q) are dom- 
inated by rare (i.e. having exponentially vanishing probability) sam- 
ples, causing Pn,j{.q) to be different from Pj^(q). We will discuss this 
issue with our numerical simulations in the next section. 

Let us start with the case T <T C . In reference [H] the computation 
of F(q) has been performed using the method of coupled replicas in the 
glassy phase (T < T c ) close to T c , with the result that, in leading order 
in both r = (T c 2 - T 2 )/2T C and q - q E A, on has F = l/6(g - q E Af- 
Let us now perform the computation of the prefactor of the cubic term 
to the next order in r. In order to compute F(q) one has to replicate 
n times both the spins o"j and the spins Tj. Two order parameter 
matrices appear: 

Qab = {o-aO-b) = (r a n) and P ab = {cr a T b ) , a,b = l, ...,n . (3) 

One makes a Parisi ansatz [2] for both Q and P, introducing two 
functions q(x) and p(x). We now recast the two nx n matrices into a 
single 2n x 2n bold-faced matrix 



Q= It r n ■ W 




Let us consider in all generality any system that is described in the 
un-coupled case by the free energy functional F[Q] in terms of the 
usual replica matrix. The same system when coupled according to J2J) 
has a free energy functional of the form 

n 

F 2 [Q,P]=F[Q]+e(nq-Y / Paa) , (5) 

o=l 

where e is a Lagrange multiplier associated to the 5- function in (j2J. 

If the problem of a single, un-coupled, system admits a continuous 
solution qp{x) (that we suppose to be known), it can be shown that 
the corresponding variational equations related to the coupled problem 
admit two simple solutions 0. 
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A first solution has AF = e = 0, with functions q(x) and p(x) that 
can be constructed explicitly from qp{x): 



where the point x is defined by the equation qp( x ) = Q- 

In the second solution q(x) = p(x). This second solution becomes 
degenerate with the first one at q = qEA- If we take the usual replica 
approach attitude to maximize the free energy with respect to the Q- 
parameters (in agreement with the physical intuition about the prob- 
lem) we will select the first solution in the region with q < qEA and 
the second one in the region with q > qEA, where it has a larger free 
energy. 

Notice that while the first solution has a flat free energy F(q) = 0, 
the second, as a consequence of the facts that 

• it coincides with the first one for q = qEA] 

• it has to be a solution of the saddle point equations; 

must be such that dF J^ \ q=qEA = 0. All these properties together 
imply that generically F(q) = A(q — qEA) 3 , where A is a problem and 
temperature dependent constant. Of course the accidental vanishing 
of A is an allowed possibility, and in that case F(q) would be of order 
five or higher. 

Let us now turn to the explicit computation for the SK model close 
to T c . In this case the replica free energy density can be written as 



F 2 [Q,P]=r TrQ* 2 + -TrQ* 3 + -y £ Q^ 4 + e(nq - £ P* aa ) , (7) 



where y = 2/3, but will be kept as a parameter during the com- 
putation, and we have defined Q* b = Q a .b/T 2 ^ QabO- + 2r) and 
p ab = p ab/T 2 . Here the solution can be found explicitly, writing 
p* = P* a and inserting P* b = Q* ab (a / b) in the equation of motion 
one finds that the matrix Q* verifies the equation: 



which is identical to the solution of the free problem with r — > (r + p*)/2 
and y — > y/2. Plugging then the known solution of the free problem 




(6) 



l,2n 



a=l 



(r+p*)Q: b + (Q* 2 ) ab + ^Q:i = o, 



(8) 
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into the free energy functional it is easy to derive the result: 



F(M) = Q - \ry) {q* ~ q* EA f- (9) 
Or, using the original definition of q, and setting y = 2/3, 

F(q) = Q + \r) (q - q E A? . (10) 

Let us now discuss the paramagnetic phase, where T > T c . Here 
qEA = and we have to select the solution that has always p{x) = q(x) 
for q > and p(x) = —q{x) for q < 0. Let us consider the replica 
symmetric free energy, corresponding to p = q, q(x) = sgn(q),p(x) = 
qo. For simplicity we will write as before equations valid for q > 0, 
keeping in mind that the free energy F(q) has to satisfy F(—q) = F(q). 

Defining Dy = e~? v /\/2irdy, standard manipulations lead to the 
expression 

F(q) = -^(l+q 2 + 2q 2 -2q -2qq )-Tlog2 + Tvq 

T j Dy log [e"" + e» cosh(2/3^y)] , (H) 

that has to be extremized with respect to qo and v (the Lagrange 
multiplier associated to q). This leads to the saddle point equations 

e- + e-cosh(2/y^) 
U e~» + e^ cosh(2/3^y) ' 1 ' 

\ + q 1 f D ye v smh(2(3^y) 

% 2 ~2f3^yJ V- + e-cosh(2/V^y)- 1 ' 

These equations have always the trivial solution q$ = 0, v = tanh -1 q, 
which is correct at high enough temperatures and small values of q. 
A small qo expansion reveals that a solution with go ^ is possible 
for q > T — 1 = q c (T). Notice that 

• for T < 1 the solution is non-trivial for all values of q, 

• for T > 2 the solution is always go = 0, 

• for 1 < T < 2 there is a phase transition from the solution qo = 
for q < T — 1 to a solution go 7^ for q > T — 1. 
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0.2 0.4 0.6 0.8 1 

Figure 1: The induced overlap q as a function of the mutual overlap q that 
we force on the system for T = 1.0, 1.25, 1.5, 1.75 and 2.00 (from right to left, 
with q = Vg for T = 1). 

The correct solution in this last domain should break the replica sym- 
metry. Close to q c (T) a small q expansion as in the previous section 
can be used even for large values of T c — T. For q c (T) — q ~ 1 we only 
consider the replica symmetric approximation. It is clear that if we 
impose q = 1 the system behaves as a single system with temperature 
T/2, and we understand therefore why q c {2) = 1. 
The (/o = solution gives a free energy 

F(g) = -^(l + g 2 )-2Tlog2 + T^-Tlog(cosh(^)) . (14) 

In figure Q we show the solution for qo as a function of q. For q = 1 we 
find that q$ takes, as it should, the value qns{T/2), i.e. the Edwards- 
Anderson value at temperature T/2 in the RS approximation. The 
functional form of qo is well approximated for all values of T by qo(q) = 
QR,s(T/2)(q — q c (T)), although small quadratic deviation from this 
form are observable. In figure ((2j) we also plot the related free energy. 
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Figure 2: The free energy —j3F(q) for the same values of the temperature 
as in figure ^ (from top to bottom). In the insert we plot, for the same 
set of temperatures again, the difference — (3(F (q) — F(q)) , where F (q) is 
computed with the trivial solution q = (from top to bottom. The difference 
is zero for T = 2 and very small for T = 1.75). 

3 Numerical versus Analytic Results 

In the following we will use our numerical results obtained from sim- 
ulations of the Jij = ±1 SK model toward two different goals, that 
we have already partially discussed in the introduction. First we want 
to have a numerical check of our analytical findings: it is clear that 
we are dealing with complex set of equations, we are working in some 
asymptotic regimes, and this implies that a numerical check is quite 
welcome. The second issue is maybe less direct but also very crucial. 
It is important to understand how well numerical simulations done on 
finite lattice encompass the infinite volume physics, and, at the same 
time, how good a control we can have on the temperature range that 
we explore. The issue of the infinite volume limit is indeed of very 
large importance, and it is crucial to verify that exotic phenomena 
(like the large deviation regime we are discussing here) are already 
well quantified in the region of lengths we can study with nowadays 
computer facilities. 

Our large scale numerical simulations are based on the numerical 
optimized Monte Carlo technique of parallel tempering ^0] • They have 
already been used in our recent paper on the problem of chaos (but 
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Figure 3: 1/N log \P^{q) iV 3^ versus g for T = 0.4 : annealed (F(q), 

represented by lines) and quenched {F(q), represented by points with error 
bars) estimates. The vertical arrow indicates the value of qEA at T = 0.4. 



for the addition of the data obtained on a N = 128 site lattice). 

We have analyzed 1024 different disorder realizations for N = 
64, 256 and 1024 sites, 256 for N = 4096 and 8192 for N = 128. 
We go down to a minimal temperature value T m = 0.4 T c = 0.4, 
and we use 0.4 10 6 iterations for thermalizing and 10 6 iterations for 
measurements. We normalize probability distributions according to 
j\P N (q)dq = l. 

The first question we discuss is whether it is appropriate to measure 
the quenched average 

0F{q) = N~ l log P N ,j(q), (15) 

(that is the quantity computed using the two-replica method) or if we 
can be happy by measuring the annealed average 

0F(q) =N~ 1 logP^J{q), (16) 

(a quantity that is much easier to measure with Monte Carlo simu- 
lation). In section [21 we have discussed the issue, arguing that these 
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two averages could differ because of rare samples. We now turn to our 
numerical data in order to answer this question. 

In figure El we show results for both the quenched (equation II 5|) 
and the annealed ( equation I16|) averages at T = 0.4 as a function of 
q. We plot data for N = 64, 128, 256 and 1024. In order to make the 
figure readable: 

1. we have subtracted from the different measurements a term 1 /3 log N, 
that separates the points from different system sizes; 

2. we have not drawn the F(q) data as points with statistical error 
bars, but we have drawn a line through the data points. 

The (few) data points for F(q) are drawn with statistical error bars 
(estimated from the variance of sample to sample fluctuations). It 
appears that the data for F(q) are confined to a much smaller q range 
than the data for F(q). This is expected since for a given value of q, 
F(q) is well defined only if (the estimate of) Pjv(?) is larger than zero 
for every sample, whereas F(q) is well defined as soon as the average 
-P/v(<?) is larger than zero. 

FigureEJshows that for large system sizes the annealed and quench- 
ed estimates are consistent. The same conclusion holds for all the 
other values of the temperature that we have considered. Because of 
that we will only consider in what follows the annealed averages -P/v(g) 
and F(q). We will estimate statistical errors on our observables by a 
jack-knife procedure |llj . 

Our second (and main) numerical goal is to verify the predic- 
tions of equation EH that predict asymptotically a cubic behavior of 
l/N log Pn(i) in Q — Qea, with a coefficient computed at order r. Now, 
is all that correct? Can we already observe this behavior on the lattice 
size that we are able to thermalize in the spin glass phase? 

We report our results at different values of the temperature (T = 
1.0, 0.8, 0.6 and 0.4) in figures H El El and respectively. With the 
dotted line we plot the theoretical prediction of equation ^1 With 
an horizontal arrow we give the analytic value of PN(q — 1) (see the 
discussion later on). 

All figures show that results for the systems of different sizes are 
converging to a common (iV independent) function: our results for the 
excess free energy are indeed in the asymptotic regime. The careful 
reader will notice however that only the limit lim^oo l/N log PN(q) 
is predicted. There is accordingly an order l/N ambiguity in the 
vertical offsets of our plots for l/N log P/v(</)- Plotting for example 
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-(1/6+T/2) (x-q EA ) 
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Figure 4: 1/N log Pn(q) versus q for T — 1: points with error bars are numer- 
ical data from different lattice sizes, the dotted line is from the perturbative 
estimate in the spin glass phase. An horizontal arrow indicates the exact 
q = 1 value. The full line (labeled "RS") is from equation fTTl 
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Figure 5: As in figure 0] but with T = 0.8. The additional vertical arrow 
indicates the value of q^A- 
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Figure 6: As in figure 03 but with T = 0.6. The additional lower dotted curve 
is drawn using a coefficient modified by hand, accounting for the renormal- 
ization of the temperature when r grows. 
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Figure 7: As in figure El but with T = 0.4. 
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Figure 8: As in figure El but with T = 1.3. 

1/JVlog (Piv(g)iV- 1 / 3 ) would result in a slightly different looking fig- 
ures. 

At T = 1, where qEA = 0, the agreement of our data with the 
theoretical prediction is very good at least up to q = 0.6. Also at 
T = 0.8, where qEA is slightly larger than 0.2 (i.e. still very small) 
the agreement is very good in a very large q range (up to q ~ 0.8). In 
all these analysis we have used the value of qEA computed by Crisanti 
and Rizzo by numerical integration |12j of the Parisi solution of the 
SK model (e.g. we are not using qEA computed to order r). 

At T = 0.6, where qEA is of order 0.5 and a small r expansion with 
two terms starts to be inappropriate, the agreement starts to be less 
good. To take into account the fact that here r is large and there is 
an effective renormalization (when T = 0.6 the first order correction 
in equation 1101 is equal to 0.16 while the leading term is equal to 1/6) 
we have added a further curve where we substitute the prefactor we 
have computed analytically in equation with some effective value, 
that optimizes the matching of the analytic form to the numerical 
data: this improves very much the agreement. This effect is even 
more dramatic at T = 0.4, where again a cubic dependence with 
a renormalized prefactor fits very well the data in a large q > qEA 
range. 
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Our third numerical goal concerns the value of of Pn(q) at g = 1 
(that we plot with an horizontal arrow in the four figures). We use 
the simple relation 

llogP JV ( g =l) = -|(^T-Jr) , (17) 

between Pn(q = 1) at temperature T and the free energy of the SK 
model at temperature T and T/2. We take for Tt the infinite volume 
free energy of the model derived with high precision in |12j from Parisi 
solution. 

One sees from figures 141 71 that the prediction of equation El f° r 
Pn(q = 1) is smaller than the value that formula 1101 takes at q = 1: 
obviously this had to be expected, since equation EH is the first term 
in an expansion in powers of q — qEA, and is not supposed to be valid 
up to the q = 1 boundary, where higher orders cannot be neglected 
(the function is presumably singular at q = 1) . Note however that 
the interplay between finite size effects and higher orders in T c — T 
is not trivial. It could be studied numerically using the multi-overlap 
algorithm which allows to measure Pn(q) with high accuracy up 
to q = 1. However parallel tempering numerical data, when pushed to 
very high statistics, confirm already in a remarkable way the analytic 
value (this can be only done for small lattice sizes). At T = 1 we 
are able to push the N = 64 lattice up to q = 1 and to keep under 
control the limit of the N = 128 lattice. At lower T values the check 
is slightly more difficult, but one sees that the q = 1 limit of larger 
and larger lattices approaches better and better the analytic value. 

That the situation is well under control can be shown nicely for 
T > 1 by comparing the numerical data, in the whole q S [0, 1] range, 
with the results of our RS computation of equation II II This is done in 
figures for T = 1.0 and |§] for T = 1.3 (the highest temperature alas 
in our parallel tempering simulation). The first figure shows the full 
line (labeled "RS" ) cleanly overshooting the dotted line, to reach the 
q = 1 axis slightly above (as expected) the exact q = 1 limit 3 . Figure 
IH1 show also again the agreement our our numerics with the replica 
calculation. 

Equation ITUl suggests |Hj that (for T < T c ) Pn(q) scales asymp- 
totically as 

P N (q) ATS =F(N(q- q E A) 3 ) for q > q E A , (18) 

3 This strongly indicates that the correct large q behavior in figures 0] and is given by 
the N = 64 and 128 data, and that the rightmost N = 256 points are somewhat misleading 



14 




Figure 9: Scaling plot of Pn(q) as a function of (q — q m ax(N)) 3 N for 

T = 1.0, compared with the coupled replica estimate (with arbitrary normal- 
ization). Our estimate of q m ax{N) is reported in the figure under the name 
qbar (it is obviously zero at T c ) . 
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Figure 10: As in figure but for T = 0.8. Coupled replica estimate is 
plotted with and without the term linear in r. Note that q max {N) is not well 
determined on small systems at this temperature, that is close to T c . 



with J-(z) oc exp(— Az), with some positive constant A, for large z. If 
the above formula holds down to the maximum of Pn(q) this implies 
that the location of the peak of the probability distribution of the 
overlap on a lattice of size 7Y, that we denote by q max (N), has to 
behave like 

_i 

q m ax{N) = qEA + cN 3 , 

a behavior that has indeed been verified with high accuracy in former 
work 0. 

In what follows we will accordingly look for a scaling law in the 
form 

P N (q) 7V-3 =g(N(q- q m ax(N)Y) (19) 

using the numerical estimates obtained for q max (N)- Figures IH1 [Till fTTI 

and [Tulare scaling plots of log (PN(q)N~ 1 / 3 ^ versus N(q — q max (N)) 3 . 
The figures clearly exhibit the predicted cubic limiting behavior. At 
T = T c the data are furthermore in excellent agreement with the 
predicted slope. At T = 0.8 the term linear in r computed in this 
paper is essential in order to obtain good agreement with the data. 
T = 0.6 and T = 0.4 are clearly too far from T c in order to obtain 
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Figure 11: As in figure ITU1 but for T = 0.6. The additional dotted line is 
drawn using the same hand modified coefficient as in figure |H1 
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Figure 12: As in figureEJbut for T = 0.4. The additional dotted line is drawn 
using the same hand modified coefficient as in figure [7| 
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a good result for the slope with only two terms in the expansion in 
powers of r. Scaling is violated on small systems for the rightmost 
points. It must be so since those are points with q close to one, and 
we have seen before that close to one the behavior is not cubic. 

It should be stressed that with our data scaling in ((/ — Qmax 
(equation ITH|) is much better than scaling in (q — qEA^N (equation 
ITT)1 using the very precise data for qEA by Crisanti and Rizzo). One 
reason for this better agreement may be that for q > qEA the inflec- 
tion point of (q — q m ax(N)) 3 N mimics the maximum of 1/N log Pj^(q). 
Scaling plots as function of (q — qEA^N with an ad-hoc effective N 
independent qEA appear (in a small z interval) in [7j for T = 0.8, and 
in (H| for T = 0.5. 

4 Conclusions 

The results of these analysis are very positive. We have first argued 
that the large deviation leading behavior of Pn(q > Qea) is gener- 
ically 1/N log Pn(i) = — A (q — qEA) 3 , and we have computed the 
first corrections to the expansion of A in powers of T c — T for the 
Sherrington-Kirkpatrick model. 

Data from numerical simulations confirm that this behavior can 
be detected already on moderate lattice sizes. First we have seen that 
we are well allowed to take an annealed average of F. Secondly we 
have analyzed the details of the cubic behavior, and checked the valid- 
ity of the perturbative estimate of the prefactor: close to T c it works 
very well, while far away from T c a simple renormalization improves 
the agreement in a substantial way. Our very accurate data have al- 
lowed a detailed analysis of the P(q = 1) value, allowing again for a 
positive check of the analytic result. For T > T c an exact computa- 
tion in the replica symmetric scheme shows nicely the crossover from 
the q « qEA behavior to the ultimate limit q = 1. At last we have 
shown that the scaling plots for Pjy (q) improve appreciably if one uses 
(q — QmaxiN)) 3 N as the scaling variable, showing clearly the predicted 
cubic behavior of log (P/v(g) iV" 1 / 3 ) . 
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